% code for robust averaging of SPMs
% Ref:1) Ritter, O., A. and Dawes, G.J.K., 1998. New equipment and processing 
%      for magnetotelluric remote refernce observations, GJI. 132, 535-548.  
% 
% 		2) Jones et al 1989, JGR.
		
% NOT tested
% development history
% 20.3.3 coding according to Oliver Ritter, 1998.
% 5.4.3 initial tf is now computed from averaged SPM
% 22.4.3 initial tf is computed from medain SPM
% 23.4.5 Tukey weights added to reduce extreme outliers (Beaton & Tukey, 1974)

% 5/6/3 trying co comapre the effect of influence functions
% Latest date 05.06.2003


function[tf_t,ProcDef] = robspm1(SPMall,ProcDef),

MinC = 'Ex';
A = size(SPMall);
Influence = ['Hubr';'L1L1';'L1L2';'LpLp';'Fair';'Cuch';'Germ';'Wels';'Tuky';'L2L2'];
SPMall = NormSPM(SPMall); 		% Normalizing SPM wrt power in horizontal magnetic fileds

%SMed = MedSPM(SPMall); 			% median auto & cross spectra
%tfmed = tf(SMed,ProcDef); 		% tf for all frequency from median spectra

global c;
c = 1.39998;



for i = 1:A(1),
   
   SPMat(:,:,:) = SPMall(i,:,:,:); % all the events for one frequency
   tfmed = tf(SPMat,ProcDef);
   ti = median(tfmed); 						% single Z matrix 
  
for ii = 1:A(2),
	data(:,:) = SPMat(ii,:,:);  		% copy of one event all frequencies
   S(ii) = abs(r_resid(ti,data,MinC)); % residual to average tf
end;											% the residual returned is |EobsEobs* - EpreEpre*|

sig_m = 1.483*median(abs(S-median(S)));% scale

X = S/sig_m; % normalize the residual w.r.t scale
for ii = 1:10,
w_m = eval([Influence(ii,:) '(X)']);
ProcDef.Weights(i,:) = w_m;
[SPMi] = WeightSPM(SPMat,ProcDef,i); % weigthing matrices
tf_t(ii,i,1:4)  = tf(SPMi,ProcDef);
tf_t(ii,i,5:8) = 0;
end;
   
   
end;


%helper functions

function[w] = L2L2(S),
d = length(S);
w = ones([1,d]);

function[w]= L1L1(S),
w = 1./abs(S);

function[w] = L1L2(S),
w = 1./(1+S.^2/2);

function[w] = LpLp(S),
w = abs(S).^(1.2-2); % v=1.2 gives a good balance between L1 (v=1) and L2(v=2)

function[w]= Fair(S),
global c;
w=1./(1+abs(S)/c); %c = 1.3998 gives 95% efficiency on Gaussian data

function[w]=Cuch(S),
global c;
w=1./(1+(S/c).^2);

function[w] = Germ(S),
w = 1./(1+S.^2).^2;

function[w] = Wels(S),
global c;
w = exp(-(S/c).^2);

function[w] = Tuky(S),
global c;
 k = find(S <= c);		% new weights
 w(k) = (1-(S(k)/c).^2).^2;
 q = find(S > c);
 w(q)=0;



 function[w]= Hubr(S),
  k = find(S <= 1.345);
 w(k)=1;
 q = find(S>1.345);
 w(q) = 1.345./abs(S(q));
 


   

   
   
   
   
   
      
      
